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Abstract 



' A lattice Boltzmann scheme able to model the hydrodynamics of phase sep- 



aration and two-phase flow is described. Thermodynamic consistency is en- 
sured by introducing a non-ideal pressure tensor directly into the collision 



^ ' operator. We also show how an external chemical potential can be used to 

(N ■ 

O ' supplement standard boundary conditions in order to investigate the effect 

o : 

' of wetting on phase separation and fluid flow in confined geometries. The 

in 

Q\ ' approach has the additional advantage of reducing many of the unphysical 

' discretisation problems common to previous lattice Boltzmann methods. 
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The hydrodynamics and kinetics of two-component fluids present a wealth of physical 
problems of both fundamental and technological importance [^. There is much current 
interest in the relevance of hydrodynamics to spinodal decomposition and the effects of 
substrates with different wetting properties on phase separation and domain growth 0. In 
addition, the flow properties of multicomponent systems, particularly in porous media, have 
been intensively studied and are of great relevance to oil recovery f^. 

Conventional methods for simulating two-phase flow include numerical integration of 
the Navier-Stokes equations and molecular dynamics simulations [^. These techniques are 
extremely computationally intensive and particularly difficult to implement in random ge- 
ometries. A newer approach, the lattice Boltzmann method, has recently proved competitive 
0. Here a set of distribution functions defined on a lattice is allowed to relax to equilibrium 
via a Boltzmann equation, discrete in both space and time. The correct choice of equilib- 
rium distribution ensures that in the long wavelength limit the Navier-Stokes equations are 
recovered. 

Several authors have set up lattice Boltzmann schemes for two-phase systems. In most 
approaches interface formation has been introduced phenomenologically by modifying the 
Boltzmann collision operator to impose phase separation |^. Recent work by Shan and 
Chen [H has attempted to relate phase separation to microscopic interactions by redefining 
the equilibrium velocity distribution so as to simulate a fluid with a non-ideal equation of 
state. However, their approach leads to inconsistent thermodynamics unless a particular 
equation of state is chosen. In addition, all current schemes reach equilibrium distributions 
which have unphysical velocity fluctuations within the interfacial region P]. 

In this letter we show for the first time that it is possible to set up a lattice Boltzmann 
scheme modelling isothermal hydrodynamics for two-phase systems. This is achieved by 
introducing directly into the collision operator the equilibrium pressure tensor for a non-ideal 
fluid. The resulting phase transition is pressure driven, as pertinent to a liquid- vapour system 
quenched to well below the critical point . The fluid reaches the correct thermodynamic 
equilibrium as determined by the equation of state and a Maxwell construction. 
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We first summarise the relevant results from the Van-der-Waals formulation of quasi- 
local thermodynamics for a two-component fluid in thermodynamic equilibrium at a fixed 
temperature [|lT|. The free energy functional is taken to be 

^(r) =y"|||Vn(r)p-hV^(n(r))|rfr, (1) 

where the first term gives the contribution from any density gradients and the second de- 
scribes the bulk free energy density. The non-local pressure is defined by 

p(r) = n— '^{r) = Po — KnW\ |Vn|^, (2) 

where po = nip'{n) — ip{n) is the equation of state of the fluid. To obtain the full pressure 



tensor in a non-uniform fluid, non-diagonal terms must be added [|T2 



P.,(r)=p(r)5., + .^:^. (3) 

In equilibrium the components of the pressure tensor define the surface tension in inhomo- 
geneous regions of the fluid. 

We next describe how a lattice Boltzmann simulation can be set up to ensure that these 
equilibrium conditions are met while incorporating the dynamical behaviour pertinent to 
fluids, namely, order-parameter conservation and hydrodynamic transport. To illustrate our 
method, we choose to work in two dimensions on a triangular lattice, the simplest geometry 
that allows us to reproduce the Navier-Stokes equations. Let fi{x,t) be a non-negative real 
number describing the distribution function of the fluid density at site x at time t moving 
in direction Cj, i = 1 . . . 6. The unit vectors Cj = {cos(27r(z — l)/6), sin(27r(z — l)/6)} are the 
velocity vectors along the links of the lattice. With each site is also associated a function 
fo{x,t) which corresponds to the component of the distribution with zero velocity. 

The distribution functions evolve according to a Boltzmann equation which is discrete 
in both space and time 

fi{x + e^, t + 1) - fi{x, t) = Qiix, t). (4) 
The most convenient choice for Qi is a single relaxation time form 



^i---{f^-fn■ (5) 
T 

The density n and macroscopic velocity u are defined by 

n^T.^ (6) 

i 

nUa = Y^fieia-, (7) 

i 

and the equihbrium distribution, /j^^, is chosen so as to reproduce the correct dynamic 
equations for n and u. 

For a one-component fiuid with an ideal gas equation of state, /^^^ is expanded as a 
power series in the local velocity and the coefficients determined by local conservation of 
mass and momentum and by the constraints of Galilean invariance and isotropy of the 
pressure tensor. In order to include the correct non-local thermodynamic properties of a 
non-ideal fiuid, additional non-local terms are needed in the expansion for f^'^. We thus 
define 

f^'i = A + BCiaUa + Cm^ + DUaUf^eiaCifS + F^Cia 

+G (8) 

f,'^ = A, + Cou\ (9) 

where the coefficients, now functions of n and its derivatives, can be determined by three 
macroscopic constraints. 

The first two are, as for the non-interacting case, local conservation of mass and momen- 
tum 

E/r=^, (10) 

i 

Yl fi'^^ioc = nUa- (11) 
i 

The third constraint is that the pressure tensor takes the form 

J2 fr^ia^^m = PccP + nUaUjj. (12) 



The constraints (|TOD, ([TlD and ([T^), together with the equihbrium thermodynamic defi- 
nitions (l)-(3), are sufficient to determine the coefficients in the expansions (||) and (|^) 

Aq = n — 2(po — KnW'^n) 
A = {pq- /tnV^n)/3 

B = n/3 Co = —n 

C = -n/6 D = 2n/3 



K 



G^y = 2-1^1^ (13) 
^ 3dxdy ^ ^ 

where, on the lattice, derivatives are simply expressed by finite difference approximations. 

The continuum hydrodynamic equations modelled by this dynamic scheme can be de- 
termined by performing a Chapman- Enskog expansion on the Boltzmann equation (1). To 
second order, the usual continuity and Navier-Stokes equations result 

^ + ^ = 0, (14) 

at OXa 

+ = -8^^+''^ + d^. (AH V . nn) , (15) 



where u = (2r — l)/8 and 



Note that to this order, the only difference between equation ([15|) and the Navier-Stokes 
equation for an ideal fluid is the appearance of the non-ideal pressure po. 

To test the correctness and applicability of the approach described above we performed 
simulations on a Van-der-Waals fluid for which 



%jj = nT In ( ) — an^. (17) 



We first consider the equilibrium configuration for a system with periodic boundary con- 
ditions and zero net flow velocity. The inset of Figure 1 shows the coexistence curve as 
a function of temperature T, calculated from equation (2) using a Maxwell construction. 
The points show simulation data obtained from lattices of size 256 x 256, equilibrated for 
10, 000 time steps. Because the correct equilibrium thermodynamics is inherent in the model, 
the bulk phases reached in the simulations obey the Maxwell construction. Note the wide 
range of coexisting densities that can be reached by the simulation before significant finite 
difference errors come into play as the temperature is lowered. 

Figure 1 also shows the equilibrium interface density profiles which are seen to depend, 
as expected, on k and the temperature. Again the agreement with a direct integration of 
the continuum thermodynamic equations is excellent. The interfacial width can be varied, 
typically between ~ 2 — 30 lattice sites. This ensures that lattice anisotropy effects can 
be made unimportant. Our mechanism for interface formation also reduces the magnitude 
of the microscopic velocities in the interfacial region by a factor of ~ 10'^ relative to other 
non-local models. This is because in equilibrium, the effective pressure force goes to zero 
and thus does not need to be balanced by momentum transfer on the scale of the lattice 
discretisation. 

In Figure 2 we emphasise the consistency between the mechanical definition of surface 
tension 

a = RAP, (18) 

where AP is the pressure difference between the inside and outside of a spherical domain of 
radius R, and the thermodynamic definition 



. = (19) 

for a flat interface. Agreement is seen as i? ^ oo with the expected curvature correction to 



equation ([T8|) appearing for i? ~ 10 lattice units. 

Far from the interface, the fluid obeys the usual Navier-Stokes equation common to other 
lattice Boltzmann schemes. However, the dynamical behaviour of the interface itself is of 



importance if domain growth is to be correctly described by the model. In Figure 3 the 
dispersion relation for capillary waves is displayed |T^ giving a best fit of u; ~ k^'^. We 
attribute the slight discrepancy from the expected dispersion relation uj"^ ~ to curvature 
corrections to Laplace's Law. These results were obtained by imposing a sine-wave of given 
wavevector on an interface that had reached equilibrium in a 128 x 128 system and observing 
the period of the subsequent oscillations for, typically, 500 timesteps. 

Finally, we demonstrate how the addition of an external chemical potential at the surfaces 
of a confined system can be used to supplement the usual bounce-back boundary conditions, 
allowing us to change the substrate properties and hence study wetting. Gradients in the 
chemical potential jj^exiz) act as a thermodynamic force on the fiuid and can be included 
within the lattice Boltzmann framework by modifying equation(ll) 

Yl fi''^i» = - ttS^ (20) 
which, in turn, introduces a force into the collision operator, equation(8), 

Fa = (21) 

As an illustrative example, by letting ^exiz) differ from zero only at the boundary sites, 
the affinity of the boundaries for each of the phases can be tuned in a simple, physically 
appealing way. Results in Figure 4 show how the fiuid configuration develops in a one- 
dimensional pore for a situation when the black, dense, fiuid (a) wets and (b) does not wet 
the walls. By modifying the functional form of Hexix), for example by introducing long range 
interactions, equilibrium phases comprising of plugs, tubes and capsules discussed by Liu 
et. al. [|l^] can be reached. In addition, an important consequence of this formalism is that 
the fiuid-boundary interface is diffuse. This allows for a reduction of many of the problems 
associated with surface orientation and standard bounce-back boundary conditions |16| . 



To summarise, we have described a lattice Boltzmann scheme, the main new features of 
which are the direct introduction of a non-ideal pressure tensor and an external chemical 
potential. This enables us to obtain an isothermal model of phase separation which correctly 



describes bulk and interfacial dynamics at low temperatures. The method also provides a 
convenient, physically motivated way of tuning boundary conditions, giving a new approach 
to situations when flow and phase separation are affected by fluid-substrate interactions. 
Moreover, unphysical velocity oscillations at surfaces and interfaces are substantially re- 
duced. 

The simplicity of the method and the ease of implementation suggests that our approach 
may be a valuable tool in the study of multi-phase hydrodynamical systems. Furthermore, 
the introduction of a temperature parameter in our mechanism for interface formation should 
allow the formalism to be extended to non-isothermal situations where heat transfer is 



important [17]. We are currently in the process of investigating such a model. 



We are indebted to Peter Coveney, Bruce Boghosian and Tim Newman for many stimu- 
lating discussions. 
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FIGURE CAPTIONS 



Fig. 1 Equilibrium density profiles normal to a flat interface for a Van der Waals fluid for 
the three highest values of T shown on the coexistence curve (inset). The solid lines 
are numerical solutions of the continuum thermodynamic equations while the points 
are from the lattice Boltzmann simulations. The parameter values are a = 9/49 and 
b = 2/21, while k, = 0.01 for the bold curves and 0.02 for the dashed curve. 

Fig. 2 AP plotted vs. 1/i? as a test of Laplace's Law. The solid line has a gradient 



calculated using equation (p^) 



Fig. 3 The dispersion curve for capillary waves on an interface, plotted on a log-log scale. 
The best fit line has gradient 1.6 ± 0.05. 

Fig. 4 Time evolution (vertically downwards) of phase separation in a narrow capillary. In 

(a) the dense (dark) fluid wets the surfaces while in (b) non-wetting is illustrated. The 
simulations were performed on a 128 x 32 lattice at T = 0.56 with (a) Hex = —0.1 and 

(b) fJ^ex = 0.5. 
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